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Abstract 

An algorithm is presented for the best least-squares 
fitting correlation matrix approximating a given missing 
▼alue or ioproper correlation matrix. The proposed algorithm 
is based upon a solution for Hosier's oblique Procrustes 
rotation problem offered by Ten Berge and Nevels. It is shown 
that the minimization problem belongs to a certain class of 
convex programs in optimization theory. A necessary and 
sufficient condition for a solution to yield the \mique 
global minimum of the least-squares function is derived from 
a theorem by Shapiro. Empirical verification of the condition 
indicates that the occurrence of non-optimal solutions with 
the proposed algorithm is very unlikely. 

ley words: missing value correlation, tetrachoric correla- 
tion, indefinite correlation matrix, constrained 
least-squares approximation, semi-infinite pro- 
gram, convex program. 
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Laast-squares Approximation of an Improper a Proper 
Correlation Matrix Using a Semi-infinite Program 

When product-fiooment correlations of a set of n variables 
are computed by any of the missing value correlation methods 
described by Frane (1978). it may happen that the resulting 
missing value correlation satrix is indefinite, and hence 
lOBpropar. This can be a serious problem in various multi-- 
variate data analysis techniques, e.g.. in regression and 
factor analysis. 

One possible approach to this problem consists of 
avoiding an (indefinite) improper correlation matrix entirely 
by estimating the missing data themselves. Hissing data can 
be estimated by maximum likelihood estimation from incomplete 
data (Beale & Little. 1975; Dempster. Laird 4 Rubin. 1977; 
Orchard & Woodbury. 1972) and by pragmatic procedures (Frane. 
1976. 1978; Gleason & Staelin. 1975; Titnn. 1970). 

Another possible approach to the problem is to render 
the improper correlation matrix non-negative definite by some 
smoothing procedure (Devlin. Gnanadesikan & Xettenring. 1975. 
p. 543: Dong. 1985; Frane. 1978). 

The purpose of the present paper is to offer a 
leasts-squares smoothing procedure. That is. one may seek the 
best fitting (in the sense of least-squares) symmetric, unit- 
diagonal, non-negative definite matrix G to the given 
improper missing value correlation matrix R. Specifically, 
the function 
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(1) 9(6) • « tr (G - R)2 

can bo minimized siibject to the constraints G = G' , 
Diag (G) = 1^ and G i 0. For convenience we wri*;e Y i 0 and 
T > 0 to denote that a symmetric matrix T is non-negative 
definite and positive definite, respectively. 

The minimization problem (l) can be generalized in three 
vays. Firstly, the problem can be applied to any improper 
correlation matrix, e.^., an indefinite tetrachoric corre- 
lation matrix or a correlation matrix obtained by element- 
wise robust estimation (Devlin, Gnanadesikan & Kettenring, 
1975. 1931. Gnanadesikan & Kettenring, 1972). Secondly, the 
problem can be generalized to handle indefinite matrices with 
fixed diagonal elements not necessary equal to one. For 
example, the scope of the problem can be extended to missing 
value covariance matrices with known variamces or to product- 
moment correlation matrices with known communal i ties. 
Thirdly, it is possible to exclude those product-moment 
correlations or covariances which are computed between com- 
plete variables (no missing values) from the minimization 
procedure. That is, the excluded elements of R can be held 
constant in (1). Without loss of generality these elements 
can be collected in the nj x nj (0 5 nj < n) submatrix 
K^x ^ 0 of S, where R is partitioned as 
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Rl2 



R = 



R2I 



R22 



In order to incorporate these three generalizations, we 
shall adress the generalized problem of minimizing (1) 
subject to the constraints 

(2a) G = G' . 

(2b) G 2 0 . 

(2c) Gjj = Ru i 0 



and 



(2d) 



Diag (G22) = Diag (R22) 0 . 



where G is partitioned as 



Gil 



G12 



G2I 



G22 
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and Gji is of order aj x . Note that the constraints (2c) 
and (2d) for the problems with n^ « 0 and - 1 are 
•quivalent. In the next section a computational solution will 
be offered for the generalized problem of minimizing (1) 
subject to the constraints (2). 

An algorithm 

The constraints G » G* (2a) and 6 ^ C (2b) can 
equivalently be expressed by the constraint 

(3) G » AA' 

for some n x m (n^ ^ m ^ n) matrix A. Considar the 
partitioning 



A = 







All 


Al2 


42 




A21 


A22 











where is of order z m, kn is of order z n^. and 
is fized in advance as 

(4) Ai = t Rii I 0 ] . 
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This choice of satisfies ths constraint G^i » Rn (2c) and 
can be adopted without loss of generality, because every 
matrix A satisfying (3) is determined up to an orthogonal 
rotation. 

Upon substitution of (3; and (4) for G in (1), the 
problem of minimizing (1) subject to the constraints (2) cap 
be reduced to the problem of minimizing the function 

(5) f(&2) - H tr (A2A2 - ^22)^ 

+ tr (AjAi - Ri2J'<*1^2 - R12J 

subject to the constraint Dlag (A2A2) - Dlag (R22)- 

In order to simplify the notation, let for any positive 
integer t the index set be defined by the Cartesian 
product 

u2 . (1 t) X il 0 , 

2 

and let t be the symmetric subset of defined by 
T«{(i,3):ii*aat(i,3)€uJ- uj^) . 

Then the ffllnloiization problem (5) can be written as 
minimizing 

(6) f(A ) s K Z Z (a'a - r )^ 

2 (i.3)«T i 3 ij 
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subject to th« constraints aj^'ajj = rj^j^ (k s nj+l n). 

where R ■ tr^j] and a^' is row i (i = 1 n) of A. For each 

k (k « n^-t-l n), (6) can be irritten as 



(7) £(A^) = X Z (a'a - r )2 
2 (i.k)€T i k ik 



+ X E (a'a - r )2 
(k.J)€T k J kj 



♦ X E E (a'a - r )2 



= E (a'a - r )2 + L 
(i.k)€T i k ik k 



= E (a'a - r )2 + L 
iitk i k ik k 



, (o) (o) . (o) (o) 

k k k k k k k 



f (a ) ♦ L . 
k k k 



where Lj^ is a constant with respect to ajj, Ajc^°J is the 
matrix A with row k replaced by zeroes, and rjj^°^ is column k 
of (R - Diag (R)] . 

In the context of Hosier's (1939) oblique Procrustes 
problem. Ten Berge and Kevels (1977) have -given a solution 
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for the global mlnlmuffl of f^Caj^) subject to the constraint 
*k'^k ^ 1* V^^^ minor adjustments, their solution can be 
generalized to oinimize £]c(a]c) subject to any arbitrary 
constraint aj^'a]^ » rj^ ^ 0. After taking a suitabJ^ initial 
choice for A2. and rov-wise rinimization 01 '7) for 

k s n^^-l n with the adjusted Ten Berge and Nevels 

solution, an algorithm for solving (5) is obtained. For each 

k (k s nj+l n). f(A2) decreases with the row-wise 

minimization, affecting only elements of row k and column k 
of AA* . The n2 ■ n — nj minimization 8tep5 can be repeated 
\mtil no significant decrease of f(A2) between two succeeding 
iteration cycles occurs. Because f(A' decreases 
monotonically and f(A2) is bounded below, convergence of the 
algorithm is guaranteed. In the next section we shall 
describe a necessary and sufficient condition for a global 
minimum of f v. A2} « 



A necessary and sufficient gondi tion for a global minimu m 

After minimizing f]c(a]^) with the adjusted Ten Berge and 
Nevels algorithm, there exists a Lagrange multiplier Oj^ such 
that 
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(Mulaik. 1972. p. SOS). Th* Lagrange multlpllez 8]^ can be 
evaluated directly from the equations (11), (12) and (13) in 
Ten Berge and Nevels (1977, p. ^95) for their cases 1, 2 and 
3 respectively. Kewriting (8) yields 

(»<<>' •»<<» ♦ .^ai)., - (A^o't'o' . .,r^, - = 0 
and hence 

(9) A'laj^ - A'rj^ - Oj^a^^ « 0 . 

where rj^ is column k of R. It should be noted that during the 
iteration process, (9) holds for the index k only immediately 
after the minimization of row k - n^ of A2. However, after 
convergence of the proposed algorithm, (9) holds simul- 
taneously for all k (k = nj+l n). Denote for convenience 

a so lut ion of the proposed algoritm by A . Then the n2 
equations (9) can be collected in the matrix equation 

(10) A'AA^ - A'Rg - A^e^2 * ° • 

where ^2 ■ 1^211^223 and 622 ■ Diag 9^). Trans- 

]>osing (10) and rewriting yields the first-order 2::2c?ssary 
conditions for a minimum of (5) 
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and 

(lib) (12^^-^22-^2^*22 = ° • 

It should be noted thrt the first-order necessary 
conditions (11) for a minimum of (5) have been obtained from 
standard partial differentiation of a constrained function 
(of. Luenberger. 1984. chap. 10). Additional results can be 
obtained from a reformulation of the problem In terms of a 
seml-lnflnlte convex program (Shapiro, 1985). This will be 
pursued next. 

Let n(T) denote the set of symmetric n x n matrices 
^ ■ I^ljl •^^^ ^^^^ ^Ij = 0 whenever (l.j) c t. Then the 
matrix G can be written as 

(12) G = C + X . 

where Z e 0(t) and C m Ic^^] such _ Jt c^j = 0 whenever 

(l.J) € T and cj^j = rj^j otherwise. In (12) G Is decomposed as 

the sum of a matrix C containing the (n^)^ n2 known (fixed) 

elements of G, and a matrix X containing the unknown (free) 

elements of G. Inserting (12) In (1) leads to the restatement 
of the minimization problem 

(13) g(X) ■ e(G) « H tr (C ♦ X - R)2 

subject to the constraints X e n(T) and (C X) ^ 0. 
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Replacing the constraint (C <«- Z) ^ 0 by the equivalent 
constraint 

(14) h(Z.u) - u'(C <f Z)u ^ 0 

for all u € *P • {u € H**: u'u s 1) makes problem (13) a 
semi— infinite program. 

Assuming that we have Rjj > 0 it can be verified that 
the semi-infinite program defined by (13) and (14) has the 
following nine properties: 

(PI) n(T) is convex. 
(P2) g(X) is convex. 
(P3) h(.,u) is concave for all u e 

(P4) The Slater (1950) condition (cf. Stoer & Wltzgall. 
1970. p. 247) holds, i.e. there exists a matrix 
Z € n(T). viz.. Zq b 0. such that h(ZQ.u) > 0 for all 
u € S'. 

(P5) 4* is compact. 

(P6) g(Z) is continuously differentiaDle. 

(P7) h(..u) is continuously different iable for all u € H'. 

(PB) h(Z.u) is continuous. 

(P9) gradj h(Z.u) is continuous. 

Properties (PI) through (P3) make the program a convex 
program znd properties (P4) through (P9) are regularity 
conditions. 
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For 8«ml<-ln£inlte programs satisfying the conditions 
(PI) through (P9), Theorem 2.2 of Shapiro (1985) is 
applicable, which states: A feasible € n(T), i.e. (C -t- Z*) 
^ 0, is a solution of the minimization problem if and only if 
there exists an n z n matrix B « (b^j] satisfying 

(i) B a B' . 

(11) (C + X*)B « 0 , 

(Hi) grad g(X)|xsX» = Pt^B) , 

^here Pr(B) is the projection of B onto the space 
n(T) defined by 

bi^ whenever (1,J) € t , 
IPT(B)]ij • { ^ 

0 otherwise , 

(Iv) B ^ 0 . 

In order to assess whether these neces. ary and 
sufficient conditions are satisfied after convergence of the 
proposed algorithm, we shall use the following lemma. 

Lemma 1. For the matrix 



(IS) B m V'B22^ • 
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conditions (1) through (ill) are satisfied, and condition 
(iv) is equivalent to the condition B22 0- 

2xSiStl. Condition (i) is obviously satisfied. 
To prove condition (ii), note that 

(16) B22WA = [ 0 I B22A22 1 • 

Sewrlting (lib) as B22A22 = 0 and transposing (16) yields 



(17) AA'WB22W = 0 . 

Substituting (C ♦ X*) = AA* and (15) in (17) proves condition 



To prove condition (ill), the matrix B is written out as 



A'WB22 = 0 



and hence 



(ii). 



(B22A2iAii)'A2iAii 



-(822*21*11 >' 



B = 



-^22*21*11 



B22 
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From (11a) it follows 

(19) B22A21A7} = -(A2Ai - R21) 
Inserting (19) in (18) yields 



B = 



-(A1A2 - Ri2)A2iAn 


A1A2 - R12 


^2*1 - R21 


A2A2 - R22 - ©22 



which can be written as 



(20) B = AA' R ^ e 



where 



= c + x*-R-e . 






(A1A2 - Rl2)A2lAll 


0 








0 


©22 



From (20) it is easily shown that Pt(B) = (C + X* - R). which 
equals grad g(X)|x=x*- This proves condition (iii). 
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Regarding condition (iv). it is obvious from (15) that 
the condition B ^ 0 is equivalent to the condition B22 ^0. • 

From our Leoma 1 and Shapiro's Theorem 2.2 it is obvious 
that B22 ^ 0 is a necessary and sufficient condition for a 
feasible Z*^ € Q(t) to be a solution of the minimization 
problem (13). It should be noted that, after convergence of 
the proposed algorithm. 622 can be evaluated hence the 
condition B22 ^ 0 can be verified. Moreover, when B22 ^ 0, it 
follows immediately from the strict convexity of g(X) that X* 
is the unique solution of the minimization problem (13). 
which means that if and only if B22 ^ 0 the unique global 
minimum of e(G) subject to the constraints (2) has been 
attained for G* = (C + X*). 

In the derivation of the necessary and sufficient 
condition for a solution to yield the unique global minimum 
of •(0). it has to be assumed that R^^ > 0. In the casa of 
singular R^^ the Slater condition (P4) does not hold and 
Aii"^ does not exist, hence it cannot be verified whether the 
obtained solution yields the unique global minimum of e(G) . 
However, for singular R^. Alexander Shapiro (personal 
communication. August 11. 1986) has shown that the problem of 
minimizing (1) subject to the constraints (2) can be 
transfoarmed to a problem of (lower) dimensionality [rank 
(Rll) ♦ n2l . with a (transformed) fixed submatrix R^j > 0. 
For reasons of availability, we give the proof which is due 
to Shapiro. 

Firstly, the function e(G) can be written as 
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(21) e(G) = X tr tP(G - R)P']2 « H tr (PGP* - PRP')2 

for any orthogonal matrix P of ord«r n x n. Secondly, let us 
take P in the form 



P = 



Pll 



such that 



(22) PiiRiiPii = 



0 


- 

0 


0 


Rii 



with R^^ > 0. Then the constraints C2> become 
(23a) PGP* = PG'P* . 



(23b) PGP- = 



PiiQiiPil 


P11G12 


<32lPil 


^22 



i 0 . 
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(23c) PuGnPii = PllRllPil ^ 0 
and 

(23d) Diag (G22) - Diag (R22) ^ 0 . 

Prom (22) and (23c) it follows that the first (nj - rank 
(Rll)] diagonal elements of PGP' are zero. From this and 
(23b) it follows that the first (n^ - rank (Kn)] rows and 
columns of P6P' are zeroes. Hence the problem of minimizing 
(21) subject to the constraints (23) is reduced to a problem 
of dimensionality (rank (Rn) *»- n2] < n. 

In order to verify the necessary and sufficient 
condition B22 ^ 0 for a solution G* s &A' to yield the unique 
global minimum of e(G) subject to the constraints (2). a 
computer program has been implemented yielding the solution 
of the minimization problem with the proposed algorithm and 
evaluating the smallest eigenvalue of 822* The computer 
program was run on 100 symmetric unit-diagonal indefinite 
matrices, where n ranged from 5 to 25. n^ ranged from 0 to 
min (10, n 2) and the column order m of A was set equal to 
n. With changes in each (free) element of G between two 
succeeding iteration cycles less than 10^ as convergence 
criterion, the algorithm never took more than 10 iteration 
cycles until convergence. Computation time never exceeded 1 
minute CPU time on a VAZ8650 computer. In all cases, the 
obtained solution satisfied the condition B22 ^ 0 within 
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accuracy limits. From thasa rasults. it can be concluded that 
tha proposed algorithm tends to produce the unique globally 
optimal solution. 

In the following leooma , another important property of 
the solution is stated. 

IifiainA_2. The rank of equals n if and only if R > 0. 

Proof - Suppose first that K > 0. Then. G* = R > 0. and hence 
the rank of equals n. 

Conversely, let the rank of G* = (C + X*) equal n. Prom 
condition (ii) it follows that 8 = 0. and from condition 
(iii) that grad g(X)|x=x* = 0- From this it follows that 
Pt(B) = (C + X* - R) = (G* ^ R) = 0 and hence G* = R > 0. 

In practice it seems to be true that the rank of G* is 
alweys less than or equal to the number p of positive 
eigenvalues of R. Since computation time heavily depends upon 
the column order m of A, it is advised to take m = p. A 
further reduction of computation time can be accomplished by 
setting a suitable initial value for A. A reasonable initial 
value A^^J can be based upon an eigen-decomposition of R 

R « m* . 

where X is an orthogonal matrix of order n x n containing as 
columns the normal i8ed eigenvectors of R and A is a diagonal 
matrix containing the n eigenvalues of R. Let An be the 
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diagonal matrix of order p with diagonal alaments the p 
positive eigenvalues of R, and let Ip be the matrix of order 
n X p with columns the p corresponding (normalized) 
eigenvectors. In the case nj = 0. it is advised to take as 
initial value for A 

(24) a'^^ = (Diag (IDl'^IDiag (IpApI^)]"^!^^^! 

where T is an arbitrary orthogonal matrix of order p x p. In 
the case nj > 0 one can take T such that the upper p x p 
submatrix of A^<>) is in lower triangular form and replace the 
submatrix Ajfo) by (4). For the 100 least-squares problems 
used above « total computation time could be reduced by more 
than 50% using (24) as initial value, and the condition B22 
was again satisfied in all cases after convergence of the 
algorithm. 



A numerleal ^Tflnip]^ 

As an illustration and for reasons of possible checks, 
an indefinite 6x6 matrix R of polychoric correlations 
(smallest eigenvalue -.0626) published by De Leeuw (1983, p. 
121) has been analyzed with various values of nj (Rjj > 0 for 
nj s 4). In order to have Rjj > 0 for nj « 5 too. the fifth 
and sixth variable have been interchanged. The matrix R is 
given in Table 1. 



rr^ 

r 
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Insert Table 1 about here 



Tablr. 2 glTes the residual matrices (6* - R) for various 
values of nj. together with the values of e(G*). Because the 
constraints (2) for the problems with nj s 0 and nj = 1 are 
equivalent, the solutions are equal. In all cases, the 
solution satisfies the condition B22 ^ 0 within accuracy 
limits . 



Insert Table 2 about here 



It can be verified that the value of e(G*) increases as n^ 
increases, as is to be expected. 

Discussion 

A monotone convergent algorithm has been constructed for 
the best least*-squares non-negative definite approximation of 
an inqproper correlation or covariauce matrix, preserving the 
diagonal elements. Additionally a verifiable necessary and 
sufficient condition for a solution to yield the unique 
global minimum of the least«^ares function has been derived. 
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Moreover, this condition tends to hm satisfied in practice. 
Thus a possibly useful alternative to existing smoothing 
procedures has been found. 

HoweTer. the solution G* is singular except in the 
trivial case R > 0 (cf . Lemma 2) hence the inverse of 0* does 
not exist. When inversion of 0* is required for a particular 
subsequent multivariate analysis, one may impose the 
additional constraint to (2) that all eigenvalues of 0 are 
greater than or equal to an arbitrary positive constant 6. An 
algorithm for the latter optimization problem is in progress, 
but is beyond the scope of the present paper. 
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TABLE 1 



De Leeuw's t jget matrix R of polychorlc correlations 
with the fifth and sixth variable Interchanged 



▼ar 


1 


2 


3 


4 


5 


6 


1 


1.0000 












2 


.4770 


1.0000 










3 


.6440 


.5160 


1.0000 








4 


.47 SO 


.2330 


.5990 


1.0000 






5 


.6510 


.6S20 


.5S10 


.7410 


1 . 0000 




6 


.S260 


.7500 


.7420 


.SOOO 


.7980 


1.0000 
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TABLE 2 



Thtt values of q(G*) and the lower-trlangular parts of the 
residual matrices (G*^ - R) using De Leeuw's target matrix, 
for various values of nj (structural zeroes omitted) 



"1 


•(G*) 


▼ar 








(G' 


- R) 










1 




2 


3 




4 


5 


0.1 


.002760 


2 




.0108 


















3 




.0011 




.0015 














4 




.0125 




.0173 - 


.0017 












5 


- 


.0063 


- 


.0088 


.0009 


- 


.0101 








6 




.0178 




.0248 


.0025 




.0286 


.0144 


2 


.002884 


3 


— 


.0012 


— 


.0016 














4 




.0131 




.0180 - 


.0019 












5 




.0067 




.0092 


.0009 




.0105 








6 




0189 




0259 


.0027 




.0296 


.0151 


3 


.002888 


4 




0132 




0181 - 


0020 












5 




0067 




0092 


0010 




0105 








6 




0189 




0259 


0028 




0296 


.0151 


4 


.003515 


5 




0083 




0116 


0016 




0133 








6 




0225 




0312 


0043 




0359 


.0185 


5 


. 004062 


6 




0243 




0349 


0057 




0403 


.0245 
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